/*
GROUP 9
This C++ project consists of functions to compute the Greeks for European options.
The formulas used have been derived from the following references:
http://en.wikipedia.org/wiki/Greeks_%28finance%29
HAUG, E. G. (2007). The complete guide to Option Pricing Formulas (2nd ed.). McGraw-Hill.


This code consists of the functions for computing each Greek, and the supporting functions called are:
Cumulative normal distribution function (CDF), Standardized normal density function (PDF), and D1 & D2
derived from the Black and Scholes Equation.

Furthermore the functions use the stdcall calling convention so they can be used by VBA in Excel.
The functions that compute the Greeks are then exported using the definition file.
*/

//We include the cmath library needed for pow(), exp() etc.
#include<cmath>

//We define pi as a constant
const double pi = 4.0*atan(1.0); 


//We write a function CDF which is the inverse CDF for the standard normal distribution
//This approximation is accurate up to 6 decimals.
double __stdcall CDF(double X)
{
	const double a1 = 0.319381530, a2 = -0.356563782, a3 = 1.781477937,
		a4 = -1.821255978, a5 = 1.330274429;
	double x = 0, k = 0; //always good to initialise variables
	double N, CDF, n;
	x = fabs(X);
	k = 1 / (1 + 0.2316419*x);
	n = (1 / sqrt(2 * pi))*exp(-0.5*x*x);
	N = 1 - n*(a1*k + a2*k*k + a3*pow(k, 3) + a4*pow(k, 4) + a5*pow(k, 5));
	if (X >= 0)
		CDF = N;

	else  CDF = 1 - N;

	return CDF;
}

//PDF is the probability density function for the standard normal distribution-N(0,1)
double __stdcall PDF(double X) {
	double PDF = 0;
	PDF = exp(-(X*X) / 2) / (sqrt(2 * pi));
	return PDF;
}

//D1 is the D1 value from the Black and Scholes Equation
double __stdcall D1(double S, double E, double r, double vol, double tau, double D)
{
	double D1 = (log(S / E) + (r - D + 0.5*vol*vol)*tau) / (vol*sqrt(tau));
	return D1;
}

//D2 is the D2 value from the Black and Scholes Equation
double __stdcall D2(double S, double E, double r, double vol, double tau, double D)
{
	double D2 = D1(S, E, r, vol, tau, D) - vol*sqrt(tau);
	return D2;
}


/*We should note at this point that Delta, Theta, Rho, Crho and Charm are different for call and put options;
The other Greeks are independent of this distinction!
In each function the following input parameters are defined: the spot price S of the underlying asset,
the exercise price E, the risk-free interest rate r, the historical volatility vol, the time to expiry tau of the contract,
the annualized dividend growth rate D, and the type of the option (1 for Call and 0 for Put). 
*/ 

//CALL AND PUT PREMIUM: calculate the option value given the input parameters 
double __stdcall PREMIUM(double S, double  E, double r, double vol, double tau, double D, int CallPut)
{
	double PREMIUM = 0;
	if (CallPut == 1) { // i.e. call option
		PREMIUM = (S*exp(-D*tau)*CDF(D1(S, E, r, vol, tau, D))) - (E*exp(-r*tau)*CDF(D2(S, E, r, vol, tau, D)));
	}
	else if (CallPut == 0) { // i.e. put option
		PREMIUM = (E*exp(-r*tau)*CDF(-D2(S, E, r, vol, tau, D))) - (S*exp(-D*tau)*CDF(-D1(S, E, r, vol, tau, D)));
	}
	return PREMIUM;
}


// FIRST ORDER GREEKS

//DELTA calculates the delta of an option given i.e.the rate of change of the theoretical option
//value with respect to changes in the underlying asset's price.
double __stdcall DELTA(double S, double E, double r, double vol, double tau, double D, int CallPut)
{
	double DELTA = 0;
	if (CallPut == 1) { // i.e. call option
		DELTA = exp(-D*tau)*CDF(D1(S, E, r, vol, tau, D));
	}
	else if (CallPut == 0) { // i.e. put option
		DELTA = -exp(-D*tau)*CDF(-D1(S, E, r, vol, tau, D));
	}
	return DELTA;
}

//VEGA calculates the vega of an option given the input paramets and measures sensitivity to volatility. 
double __stdcall VEGA(double S, double E, double r, double vol, double tau, double D) {
	double VEGA = S*exp(-D*tau)*PDF(D2(S, E, r, vol, tau, D))*sqrt(tau);
	return VEGA;
}

//THETA calculates the theta of an option given the input paramets and measures sensitivity to the passage of time.
double __stdcall THETA(double S, double E, double r, double vol, double tau, double D, int CallPut) {
	double THETA = 0;
	if (CallPut == 1) { // i.e. call
		THETA = -exp(-D*tau)*(S*vol*CDF(D1(S, E, r, vol, tau, D)) / (2 * sqrt(tau))) - r*E*exp(-r*tau)*CDF(D2(S, E, r, vol, tau, D)) + D*S*exp(-D*tau)*CDF(D1(S, E, r, vol, tau, D));
	}
	else if (CallPut == 0) { // i.e. put
		THETA = -exp(-D*tau)*(S*vol*CDF(D1(S, E, r, vol, tau, D)) / (2 * sqrt(tau))) + r*E*exp(-r*tau)*CDF(-D2(S, E, r, vol, tau, D)) - D*S*exp(-D*tau)*CDF(-D1(S, E, r, vol, tau, D));
	}
	return THETA;
}

//RHO calculates the rho of an option given the input paramets and measures sensitivity to the interest rate.
double __stdcall RHO(double S, double E, double r, double vol, double tau, double D, int CallPut) {
	double RHO = 0;
	if (CallPut == 1) { // i.e. call
		RHO = E*tau*exp(-r*tau)*CDF(D2(S, E, r, vol, tau, D));
	}
	else if (CallPut == 0) { // i.e. put
		RHO = -E*tau*exp(-r*tau)*CDF(-D2(S, E, r, vol, tau, D));
	}
	return RHO;
}

//Lamda or Leverage (elasticity):measures the percentage change in option value per percentage change in the underlying price.
// also called gearing
double __stdcall LAMDA(double S, double  E, double r, double vol, double tau, double D,int CallPut)
{
	double LAMDA = 0;
	if (CallPut == 1) { // i.e. call
		LAMDA = (exp(-D*tau)*CDF(D1(S, E, r, vol, tau, D))*S) / PREMIUM(S, E, r, vol, tau, D, CallPut);
	}
	else if (CallPut == 0) { // i.e. put
		LAMDA = (exp(-D*tau)*(CDF(D1(S, E, r, vol, tau, D)) - 1) * S) / PREMIUM(S, E, r, vol, tau, D, CallPut);
	}
	return LAMDA;
}

// SECOND ORDER GREEKS

//GAMMA calculates the gamma of an option and measures the rate of change in the delta wrt changes in spot price.
double __stdcall GAMMA(double S, double E, double r, double vol, double tau, double D) {
	double GAMMA = 0;
	GAMMA = exp(-D*tau)*(PDF(D1(S, E, r, vol, tau, D)) / (S*vol*sqrt(tau)));
	return GAMMA;
}

//VANNA calculates the vanna of an option and measures sensitivity of delta wrt changes in volatility. 
double __stdcall VANNA(double S, double E, double r, double vol, double tau, double D) {
	double VANNA = 0;
	VANNA = exp(-D*tau)*PDF(D1(S, E, r, vol, tau, D))*D2(S, E, r, vol, tau, D)/vol;
	return VANNA;
}

//VOMMA calculates the vomma or Vega Convexity of an option; it measures second order sensitivity to volatility. 
double __stdcall VOMMA(double S, double E, double r, double vol, double tau, double D) {
	double VOMMA = 0;
	VOMMA = S*exp(-D*tau)*PDF(D1(S, E, r, vol, tau, D))*sqrt(tau)*((D1(S, E, r, vol, tau, D)*D2(S, E, r, vol, tau, D)) / (vol));
	return VOMMA;
}

//CHARM calculates the charm or delta decay of an option; it measures the instantaneous rate of change of Delta over the passage of time.
double __stdcall CHARM(double S, double E, double r, double vol, double tau, double D, int CallPut) {
	double CHARM = 0;
	if (CallPut == 1) { //i.e. call
		CHARM = D*exp(-D*tau)*CDF(D1(S, E, r, vol, tau, D)) - exp(-D*tau)*PDF(D1(S, E, r, vol, tau, D))*((2*(r-D)*tau-D2(S,E,r,vol,tau,D)*vol*sqrt(tau)) / (2*tau*vol*sqrt(tau)));
	}
	else if (CallPut == 0) { // i.e. put
		CHARM = -D*exp(-D*tau)*CDF(-D1(S, E, r, vol, tau, D)) - exp(-D*tau)*PDF(D1(S, E, r, vol, tau, D))*((2*(r-D)*tau-D2(S,E,r,vol,tau,D)*vol*sqrt(tau)) / (2*tau*vol*sqrt(tau)));
	}
	return CHARM;
}

//VETA calculates the veta of an option and measures the rate of change in the Vega wrt the passage of time.
double __stdcall VETA(double S, double E, double r, double vol, double tau, double D) {
	double VETA = 0;
	VETA = S*exp(-D*tau)*PDF(D1(S, E, r, vol, tau, D))*sqrt(tau)*(D + (((r - D)*D1(S, E, r, vol, tau, D)) / (vol*sqrt(tau))) - ((1 + D1(S, E, r, vol, tau, D)*D2(S, E, r, vol, tau, D)) / (2 * tau)));
	return VETA;
}



// THIRD ORDER GREEKS

//COLOR calculates the color or gamma decay of an option; it measures the rate of change of Gamma over the passage of time.
double __stdcall COLOR(double S, double E, double r, double vol, double tau, double D) {
	double COLOR = 0;
	COLOR = -exp(-D*tau)*((PDF(D1(S, E, r, vol, tau, D))) / (2 * S*tau*vol*sqrt(tau)))*(2*D*tau+1+((2*(r-D)*tau-D2(S,E,r,vol,tau,D)*vol*sqrt(tau))/(vol*sqrt(tau)))*D1(S,E,r,vol,tau,D));
	return COLOR;
}

//SPEED calculates the speed of an option and measures the rate of change in Gamma wrt changes in the underlying price.
double __stdcall SPEED(double S, double E, double r, double vol, double tau, double D) {
	double SPEED = 0;
	SPEED = -exp(-D*tau)*((PDF(D1(S, E, r, vol, tau, D))) / (S*S*vol*sqrt(tau)))*(((D1(S, E, r, vol, tau, D)) / (vol*sqrt(tau))) + 1);
	return SPEED;
}

//ULTIMA calculates the ultima of an option and measures the sensitivity of Vomma wrt change in volatility. 
double __stdcall ULTIMA(double S, double E, double r, double vol, double tau, double D) {
	double ULTIMA = 0;
	ULTIMA = ((-VEGA(S,E,r,vol,tau,D)) / (vol*vol))*(D1(S,E,r,vol,tau,D)*D2(S,E,r,vol,tau,D)*(1-D1(S,E,r,vol,tau,D)*D2(S,E,r,vol,tau,D)) + D1(S,E,r,vol,tau,D)*D1(S,E,r,vol,tau,D)+D2(S,E,r,vol,tau,D)*D2(S,E,r,vol,tau,D));
	return ULTIMA;
}

//ZOMMA calculates the zomma of an option and measures the rate of change of Gamma wrt changes in volatility
double __stdcall ZOMMA(double S, double E, double r, double vol, double tau, double D) {
	double ZOMMA = 0;
	ZOMMA = -exp(-D*tau)*((PDF(D1(S, E, r, vol, tau, D))*(D1(S, E, r, vol, tau, D)*D2(S, E, r, vol, tau, D) - 1)) / (S*vol*vol*sqrt(tau)));
	return ZOMMA;
}


// WITH RESPECT TO CHANGES IN THE STRIKE PRICE E

//DUAL_DELTA calculates the dual delta of an option and measures the exact  
// probability of an option finishing at a particular price of E.
double __stdcall DUAL_DELTA(double S, double E, double r, double vol, double tau, double D, int CallPut) {
	double DUAL_DELTA = 0;
	if (CallPut == 1) { //i.e. call
		DUAL_DELTA = -exp(-r*tau)*CDF(D2(S, E, r, vol, tau, D));
	}
	else if (CallPut == 0) { //i.e. put
		DUAL_DELTA = exp(-r*tau)*CDF(-D2(S, E, r, vol, tau, D));
	}
	return DUAL_DELTA;
}

//DUAL_GAMMA calculates the dual gamma of an option and measures how fast the Dual Delta changes wrt strike price E.
double __stdcall DUAL_GAMMA(double S, double E, double r, double vol, double tau, double D) {
	double DUAL_GAMMA = 0;
	DUAL_GAMMA = exp(-r*tau)*((PDF(D2(S, E, r, vol, tau, D))) / (E*vol*sqrt(tau)));
	return DUAL_GAMMA;
}
